In [1]:
import numpy as np
%matplotlib nbagg
import matplotlib.pyplot as plt
import scipy.io.wavfile

In [2]:
import os
print os.getcwd()


/Users/tdrivas/Documents/Chichi Python Class/public

In [3]:
#https://github.com/chichilalescu/python-for-scientific-computing/tree/master/data

r, d0 = scipy.io.wavfile.read('data/tst0.wav')
print type(d0)
print r
print d0.shape, d0.dtype

t0 = np.array(range(d0.shape[0])).astype(np.float) / r
print t0


r, d1 = scipy.io.wavfile.read('data/tst1.wav')
print r
print d1.shape, d1.dtype

t1 = np.array(range(d1.shape[0])).astype(np.float) / r
print t1


<type 'numpy.ndarray'>
44100
(66048, 2) int16
[  0.00000000e+00   2.26757370e-05   4.53514739e-05 ...,   1.49761905e+00
   1.49764172e+00   1.49766440e+00]
44100
(65536, 2) int16
[  0.00000000e+00   2.26757370e-05   4.53514739e-05 ...,   1.48600907e+00
   1.48603175e+00   1.48605442e+00]

In [4]:
ax = plt.figure().add_subplot(111)
ax.plot(t1, d1[:, 0])
ax.plot(t1, d1[:, 1])


Out[4]:
[<matplotlib.lines.Line2D at 0x7fed12044090>]

In [5]:
sfull1 = np.fft.rfft(d1, axis = 0)
sfull0 = np.fft.rfft(d0, axis = 0)
k0 = np.array(range(sfull0.shape[0])).astype(t0.dtype) / t0[-1]
k1 = np.array(range(sfull1.shape[0])).astype(t1.dtype) / t1[-1]
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
ax.plot(k0, np.abs(sfull0)[:, 0]**2, label = '0')
ax.plot(k1, np.abs(sfull1)[:, 0]**2, label = '1')
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('intensity')
ax.set_xlabel('Hz')
ax.legend(loc = 'best')


Out[5]:
<matplotlib.legend.Legend at 0x7fed11d06d10>
\begin{equation} f: \mathbb{R} \rightarrow \mathbb{R} \end{equation}

but we don't want $10^{-1}$ Hz, we only want stuff above 20Hz. or 16, if you have really good ears.

16 Hz means 1/16 seconds, and we have a sampling rate of 44100. Therefore we want 44100/16 samples at most in a signal


In [6]:
print t0.shape[0]/(44100/16.), t1.shape[0]/(44100/16.)


23.9629931973 23.7772335601

In [7]:
nsamples = int(t0.shape[0]/(44100/16.))
nsamplestmp = int(t1.shape[0]/(44100/16.))
if nsamplestmp < nsamples:
    nsamples = nsamplestmp
print nsamples


23

In [8]:
d0 = d0[:nsamples*int(44100/16.)].reshape(nsamples, int(44100/16.), 2)
d1 = d1[:nsamples*int(44100/16.)].reshape(nsamples, int(44100/16.), 2)
t0 = t0[:int(44100/16.)]
t1 = t1[:int(44100/16.)]

In [9]:
sfull1 = np.fft.rfft(d1, axis = 1)
sfull0 = np.fft.rfft(d0, axis = 1)
k0 = np.array(range(sfull0.shape[1])).astype(t0.dtype) / t0[-1]
print k0
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
s0 = np.average(np.abs(sfull0)[:, :, 0]**2, axis = 0)
s1 = np.average(np.abs(sfull1)[:, :, 0]**2, axis = 0)
ax.plot(k0, s0, label = '0')
ax.plot(k0, s1, label = '1')
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylabel('intensity')
ax.set_xlabel('Hz')
ax.legend(loc = 'best')


[  0.00000000e+00   1.60072595e+01   3.20145191e+01 ...,   2.20259891e+04
   2.20419964e+04   2.20580036e+04]
Out[9]:
<matplotlib.legend.Legend at 0x7fed11c9a350>

In [10]:
ax = plt.figure(figsize=(6, 4)).add_axes([.15, .15, .7, .7])
ax.hist(d0[:, :, 0].ravel(),
        bins = 32,
        histtype = 'step')
ax.hist(d0[:, :, 1].ravel(),
        bins = 32,
        histtype = 'step')
ax.hist(d1[:, :, 0].ravel(),
        bins = 32,
        histtype = 'step')
ax.hist(d1[:, :, 1].ravel(),
        bins = 32,
        histtype = 'step')


Out[10]:
(array([   24.,   109.,   223.,   489.,   826.,  1455.,  2226.,  2932.,
         3151.,  2964.,  3089.,  3527.,  3522.,  4000.,  3600.,  3038.,
         3056.,  3301.,  3699.,  3431.,  2692.,  1910.,  1413.,  1450.,
         1530.,  1439.,  1039.,   818.,   692.,   509.,   398.,   836.]),
 array([-30641. , -28659.5, -26678. , -24696.5, -22715. , -20733.5,
        -18752. , -16770.5, -14789. , -12807.5, -10826. ,  -8844.5,
         -6863. ,  -4881.5,  -2900. ,   -918.5,   1063. ,   3044.5,
          5026. ,   7007.5,   8989. ,  10970.5,  12952. ,  14933.5,
         16915. ,  18896.5,  20878. ,  22859.5,  24841. ,  26822.5,
         28804. ,  30785.5,  32767. ]),
 <a list of 1 Patch objects>)

In [10]: